R Notebook

#setwd("C:/Users/erussek/forage_jsp/analysis")
setwd("/Users/evanrussek/forage_jsp/analysis")
rm(list = ls())
library(tidyr)

Attaching package: ‘tidyr’

The following object is masked from ‘package:Matrix’:

    expand

library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(zoo)

Attaching package: ‘zoo’

The following objects are masked from ‘package:base’:

    as.Date, as.Date.numeric

library(ggplot2)
library(ggpubr)

Loading required package: magrittr

Attaching package: ‘magrittr’

The following object is masked from ‘package:tidyr’:

    extract

library(RcppRoll)
library(knitr)

data <- read.csv('data/run4_data.csv')
head(data)

data <- data %>% ungroup() %>% bind_cols(s_num = group_indices(., subjectID))
data <- data %>% mutate(subj = factor(s_num))

travel_keys = data %>% ungroup() %>% select(travel_key) %>% unique()
clean_subj_data <- function(subj_data, travel_keys){
  travel_key_hard = travel_keys$travel_key[1]
  travel_key_easy = travel_keys$travel_key[2]
  
  subj_data <- data %>%select(round, phase, reward_obs, reward_true, lag, exit, start_reward, n_travel_steps,
                                    travel_key, subjectID, trial_num, s_num, correct_key) %>% group_by(trial_num) %>%
    mutate(press_num = 1:n(), round = as.factor(round), 
           rep_number = case_when(trial_num < 7 ~ 1, TRUE ~ 2)) %>%
    mutate(phase = replace(phase, phase == "Harvest", "HARVEST"))

  subj_data <- subj_data %>% group_by(trial_num) %>% 
    mutate(travel_key = replace(travel_key, travel_key == "", first(travel_key[travel_key != ""]))) %>%
    mutate(travel_key_cond = case_when(travel_key == travel_key_hard   ~ "HARD",
                                       travel_key == travel_key_easy  ~ "EASY")) %>%
    filter(!is.na(travel_key_cond))
}

# clean all the data and bind
datalist <- list()
for (i in 1:10){
  subj_data <- data %>% filter(s_num == i)
  subj_data <- clean_subj_data(subj_data, travel_keys)
  subj_data <- ungroup(subj_data)
  datalist[[i]] <- subj_data
}

Unequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorUnequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorUnequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorUnequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorUnequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorUnequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorUnequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorUnequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorUnequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorUnequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vector

#cdata <- do.call(rbind,datalist)
cdata <- bind_rows(datalist)

# make a function to plot reward on each game... 
plot_subj_reward_v_press <- function(subj_data){
  
  s_num <- subj_data$s_num[1]
  
  # reward plots which show the threshold... 
  rew_plot <- ggplot(subj_data, aes(x = press_num, y = reward_obs, group = round))  +
    #geom_rect(data = t1_data, aes(xmin = press_num - 0.5, xmax = press_num + 0.5, ymin = -Inf, ymax = Inf, fill = round)) +
    geom_point(aes(color = phase)) + facet_grid(n_travel_steps ~ travel_key_cond) + theme(legend.position = "none") + ylim(0,110) + ggtitle(paste("subj: ", s_num)) + ylab('reward collected') + xlab('button press number') + theme_minimal()
  
 # plot_an <- annotate_figure(plot, top = paste('subj: ', s_num))
  plot(rew_plot)
                             
}
for (s in 1:22){
  subj_data <- cdata %>% filter(s_num == s)
  subj_data$round_num <- as.integer(as.character(subj_data$round))

  #subj_data <- subj_data %>% group_by(travel_key_cond, n_travel_steps) %>% mutate(press_num = 1:n()) %>% ungroup()
  plot_subj_reward_v_press(subj_data)  
}

subj_data$round_num 

   [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
  [41] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
  [81] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [121] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [161] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [201] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [241] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [281] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [321] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [361] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
 [401] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
 [441] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
 [481] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
 [521] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
 [561] 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [601] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [641] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [681] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
 [721] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
 [761] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
 [801] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
 [841] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
 [881] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
 [921] 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
 [961] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 6 6 6
 [ reached getOption("max.print") -- omitted 32240 entries ]

get_trial_exits <- function(thdata){
    # get harvest
  last_phase <- last(thdata$phase)
  
  # go through data.. if last phase was harvest, remove that round... 
  if (last_phase == "HARVEST"){
    last_round <- last(thdata$round)
    # get the last reward ops...
    last_reward_ob <- last(thdata$reward_obs[!is.na(thdata$reward_obs)])
    
    #if (last_reward_ob > 8){
    thdata <- thdata %>% filter(round != last_round)
    #}
  }
  
  # now select harvest data
  thdata <- thdata %>% filter(phase == "HARVEST", !is.na(reward_obs)) # find out why reward_true has so many .na
  
  ## get either last true reward observed
  return_tbl <- thdata %>% 
                  group_by(round) %>% 
                      summarise(s_num = first(s_num), last_reward = last(reward_obs),
                          trial_num = first(trial_num),
                          n_travel_steps = first(n_travel_steps),
                          travel_key_cond = first(travel_key_cond),
                          rep_number = first(rep_number)) %>% ungroup()
  return(return_tbl)
}

#sdata %>% group_by(trial_num) %>%
#  do(get_trial_exits(.))

exit_data <- cdata %>% group_by(s_num, trial_num) %>%
  do(get_trial_exits(.)) %>% ungroup() %>% mutate(subj = as.factor(s_num))

library(ggpubr)
exit_data <- exit_data %>% mutate(round_num = as.integer(as.character(round)))
for (s in 1:22){
  p <- ggplot(exit_data %>% filter(s_num == s), aes(x = round_num, y = last_reward)) + geom_point() + geom_line() +
  facet_grid(n_travel_steps ~ travel_key_cond) + theme_minimal()  + ylim(0,100) + ylab('last reward collected') + xlab('harvest number')
  
  #plot(p)
  #ggarrange(p)
  a <- annotate_figure(p, right = "# travel steps", top = "travel step difficulty", fig.lab = paste('subj: ' , s))
  plot(a)
}

geom_path: Each group consists of only one observation. Do you need to adjust the group aesthetic?

# do some aggregating over exit thresholds... 
# just get the mean for each subject for each timepoint, collapse accross rep number...
mn_exit <- exit_data %>% group_by(s_num, n_travel_steps, travel_key_cond) %>%
  summarise(rep_exit_thresh = mean(last_reward), trial_num = mean(trial_num)) %>% 
  group_by(s_num, n_travel_steps, travel_key_cond) %>%
  summarise(exit_thresh = mean(rep_exit_thresh), trial_num = mean(trial_num)) %>%
  mutate(subj = as.factor(s_num)) %>% ungroup()

Unequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vector

# now just plot this for each subject -- draw lines to show easy - hard effect
ggplot(mn_exit, aes(x = travel_key_cond, y = exit_thresh)) + geom_dotplot(binaxis = 'y', aes(fill = subj), dotsize = 1.2) + geom_line(aes(group = subj, color = subj), size = 1.2) +
  facet_grid(.~n_travel_steps) + theme_minimal() + labs(y = 'Mean Last Reward Collected', x = 'Travel Sequence', subtitle = '# Travel Presses') #theme(legend.position = "none")
ggsave('reward_press2.png')

Saving 7.29 x 4.51 in image

ggplot(mn_exit, aes(x = factor(n_travel_steps), y = exit_thresh)) +
  geom_dotplot(binaxis = 'y', aes(fill = subj), dotsize = 1.2) + 
  geom_line(aes(group = subj, color = subj), size = 1.2) +
  facet_grid(.~travel_key_cond) + 
  theme_minimal() + 
  labs(y = 'Mean Last Reward Collected', x = 'Travel Key Condition', subtitle = '# Travel Presses') 
  #theme(legend.position = "none")
ggsave('reward_press.png')

Saving 7.29 x 4.51 in image

# plot the means for each of these...
gmn_exit <- mn_exit %>% group_by(n_travel_steps, travel_key_cond) %>% summarise(gm_thresh = mean(exit_thresh), gsd_thresh = sd(exit_thresh)/sqrt(n()))

ggplot(gmn_exit, aes(x = factor(n_travel_steps), y = gm_thresh, color = travel_key_cond)) +
  geom_line(aes(group = travel_key_cond), size = 2) +
  geom_errorbar(aes(ymin = gm_thresh - gsd_thresh, ymax = gm_thresh+gsd_thresh), width = .1, size = 2) + 
geom_point(size = 4) + ylab('group mean exit threshold') + xlab('# travel steps') + labs(color = 'travel effort') + theme_minimal()
ggsave('group mean exit thresh.png')

Saving 7.29 x 4.51 in image

library(optimx)
library(lmerTest)
exit_model <- lmer(last_reward ~ n_travel_steps + travel_key_cond  + (1  + n_travel_steps + travel_key_cond |subj), data = exit_data,  control = lmerControl(optimizer ="optimx",optCtrl=list(method='nlminb')))
summary(exit_model)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: last_reward ~ n_travel_steps + travel_key_cond + (1 + n_travel_steps +      travel_key_cond | subj)
   Data: exit_data
Control: lmerControl(optimizer = "optimx", optCtrl = list(method = "nlminb"))

REML criterion at convergence: 6956.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.1827 -0.5108  0.0002  0.5046  5.3563 

Random effects:
 Groups   Name                Variance  Std.Dev. Corr       
 subj     (Intercept)         290.72330 17.0506             
          n_travel_steps        0.06055  0.2461  -0.27      
          travel_key_condHARD  68.01256  8.2470  -0.25  0.26
 Residual                      87.82008  9.3712             
Number of obs: 928, groups:  subj, 22

Fixed effects:
                    Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)         54.95291    3.71001 20.77742  14.812 1.63e-12 ***
n_travel_steps      -0.24696    0.06318 19.11207  -3.909 0.000935 ***
travel_key_condHARD -5.06293    1.89212 17.91918  -2.676 0.015465 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) n_trv_
n_trvl_stps -0.309       
trvl_k_HARD -0.256  0.206

library(broom.mixed)
td <- tidy(exit_model, conf.int = TRUE)
tdsome <- td %>% filter(term == 'n_travel_steps' | term == 'travel_key_condHARD')

tdsome <- tdsome %>% mutate(parameter = case_when(term == 'n_travel_steps' ~ 'travel step', term == 'travel_key_condHARD' ~ 'travel effort cond'))

ggplot(tdsome, aes(estimate,parameter,color = parameter)) + 
  geom_point(size = 2)  +
geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), size = 1.2) + theme(legend.position = "none") + theme_minimal() + ylab('fixed effect') + xlab('estimate (# points)') + ggtitle('linear mixed effects model - effect on exit threshold') + theme(axis.text.y = element_text(angle = 90, hjust = .5))
ggsave('exit model.png')

Saving 7.29 x 4.51 in image

# reaction times... 


## a bit more sensible... 
lag_data <- cdata %>% select(s_num, rep_number, travel_key_cond, n_travel_steps, phase, trial_num, lag, correct_key, round, trial_num) %>% mutate(log_lag = log(lag)) %>% group_by(s_num, trial_num, round, phase) %>% slice(-1) %>% ungroup() %>% mutate(subj = as.factor(s_num))

ggplot(lag_data, aes(x=lag)) + geom_histogram(binwidth = 10) + xlim(50,300) + facet_wrap( ~ s_num)

#ggplot(lag_data, aes(x=log_lag))  + geom_histogram(binwidth = .01) + facet_wrap( ~ s_num)

filt_lag <- lag_data %>% group_by(s_num) %>% filter(log_lag < (mean(log_lag) + 1.5*sd(log_lag)) , log_lag > (mean(log_lag) - 1.5*sd(log_lag))) %>% ungroup()

ggplot(filt_lag, aes(x=lag)) + geom_histogram(binwidth = 10)  + facet_wrap(~s_num) + xlim(50,350)

ggplot(filt_lag, aes(x=log_lag)) + geom_histogram(binwidth = .02) + facet_wrap(~s_num)

## with a run, plot the mean lag 
round_filt_lag <- filt_lag %>% 
  group_by(s_num, n_travel_steps, travel_key_cond, round, phase) %>%
  summarise(trial_num = first(trial_num), mean_lag = median(lag), mean_log_lag = median(log_lag)) %>% mutate(round_num = as.integer(as.character(round)))

for (s in 1:22){
  
  # plot mean lag as a function of round number
  p <- ggplot(round_filt_lag %>% filter(s_num == s), aes(x = round_num, y = mean_log_lag, color = phase)) + geom_point(size = 2) +facet_grid(n_travel_steps ~ travel_key_cond) + geom_line(size = 1.2)  + xlab('harvest/travel number') + ylab('mean log lag') + ggtitle(paste('subj: ', s)) + theme_minimal()
  plot(p)
  
}

NA
NA

library(plotrix)
cond_round_filt_lag <- round_filt_lag %>% 
  group_by(n_travel_steps, travel_key_cond, round_num, phase) %>%
  summarise(mean_log_lag = median(mean_log_lag), sd_lag = std.error(mean_log_lag))

p <- ggplot(cond_round_filt_lag, aes(x = round_num, y = mean_log_lag, color = phase)) + geom_point(size = 2) +facet_grid(n_travel_steps ~ travel_key_cond) + geom_line(size = 1.2) + geom_linerange(aes(ymin = mean_log_lag - sd_lag, ymax = mean_log_lag + sd_lag)) +
  xlab('harvest/travel number') + ylab('mean log lag') + ggtitle('group: ') +  theme_minimal()

plot(p)

cond_round_filt_lag

trial_filt_lag <- round_filt_lag %>% 
  group_by(s_num, n_travel_steps, travel_key_cond, phase) %>%
  summarise(trial_num = first(trial_num), 
            mean_lag = mean(mean_lag), 
            mean_log_lag = mean(mean_log_lag)) %>% ungroup() %>% mutate(s_num = as.factor(s_num))

ggplot(trial_filt_lag, 
       aes(x = factor(n_travel_steps), y = mean_log_lag, fill = s_num)) + 
  geom_point(aes(color = s_num), size = 2)+ 
  geom_line(aes(group = s_num, color = s_num), size = 1.2) +
  facet_grid(phase ~ travel_key_cond) + theme_minimal() + xlab('# travel steps') + ylab('log lag time')

ggplot(trial_filt_lag, 
       aes(x = travel_key_cond, y = mean_lag, fill = s_num)) + 
  geom_point(aes(color = s_num), size = 2)+ 
  geom_line(aes(group = s_num, color = s_num), size = 1.2) +
  facet_grid(phase ~ n_travel_steps) + theme_minimal()

group_lag <- trial_filt_lag %>% group_by(n_travel_steps, travel_key_cond, phase) %>% summarise(gm_lag = mean(mean_lag), gml_lag = mean(mean_log_lag)) %>% ungroup()

ggplot(group_lag, aes(x = factor(n_travel_steps), y = gml_lag, color = travel_key_cond)) +
  geom_line(aes(group = travel_key_cond), size = 2) + facet_wrap(~phase) + geom_point(size = 2) +xlab('# travel steps') +
  ylab('group mean log lag') + theme_minimal()

NA
NA
NA

harvest_lag <- round_filt_lag %>% ungroup() %>% filter(phase == "HARVEST") %>% mutate(subj = as.factor(s_num))
head(harvest_lag)

hl_model <- lmer(scale(mean_log_lag) ~  n_travel_steps + travel_key_cond + round_num  + (1 + n_travel_steps+ travel_key_cond + round_num |subj), data = harvest_lag, control = lmerControl(optimizer ="optimx",optCtrl=list(method='nlminb')))

summary(hl_model)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: scale(mean_log_lag) ~ n_travel_steps + travel_key_cond + round_num +  
    (1 + n_travel_steps + travel_key_cond + round_num | subj)
   Data: harvest_lag
Control: lmerControl(optimizer = "optimx", optCtrl = list(method = "nlminb"))

REML criterion at convergence: 1236.7

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.6966 -0.4726 -0.0215  0.4743  7.2375 

Random effects:
 Groups   Name                Variance  Std.Dev. Corr             
 subj     (Intercept)         0.7512496 0.86675                   
          n_travel_steps      0.0001744 0.01320  -0.06            
          travel_key_condHARD 0.1057507 0.32519  -0.08 -0.28      
          round_num           0.0001980 0.01407  -0.33  0.27 -0.07
 Residual                     0.1563930 0.39547                   
Number of obs: 1029, groups:  subj, 22

Fixed effects:
                     Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)         -0.067679   0.189199 21.299098  -0.358 0.724080    
n_travel_steps      -0.000230   0.003159 21.798357  -0.073 0.942619    
travel_key_condHARD  0.076090   0.074627 19.314980   1.020 0.320517    
round_num            0.023682   0.005448 15.361270   4.347 0.000546 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) n_trv_ t__HAR
n_trvl_stps -0.128              
trvl_k_HARD -0.103 -0.220       
round_num   -0.298  0.228  0.012

library(broom.mixed)
td <- tidy(hl_model, conf.int = TRUE)
tdsome <- td %>% filter(term == 'n_travel_steps' | term == 'travel_key_condHARD' | term == 'round_num')

tdsome <- tdsome %>% mutate(parameter = case_when(term == 'n_travel_steps' ~ 'travel step', term == 'travel_key_condHARD' ~ 'travel effort cond', term == 'round_num' ~ 'round number'))

ggplot(tdsome, aes(estimate,parameter,color = parameter)) + 
  geom_point(size = 2)  +
geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), size = 1.2) + theme(legend.position = "none") + theme_minimal() + ylab('fixed effect') + xlab('estimate (# points)') + ggtitle('linear mixed effects model - effect on lag') + theme(axis.text.y = element_text(angle = 90, hjust = .5))